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ABSTRACT 

We present a new algorithm for computing hyperexponen- 
tial solutions of ordinary linear differential equations with 
polynomial coefficients. The algorithm relies on interpret- 
ing formal series solutions at the singular points as analytic 
functions and evaluating them numerically at some common 
ordinary point. The numerical data is used to determine a 
small number of combinations of the formal series that may 
give rise to hyperexponential solutions. 

Categories and Subject Descriptors 

1.1.2 [Computing Methodologies]: Symbolic and Alge- 
braic Manipulation — Algorithms 

General Terms 

Algorithms 

Keywords 

Closed form solutions, D-finite equations, Effective analytic 
continuation 

1. INTRODUCTION 

We consider linear differential operators 



P = p r D T + Pr -lD r 



+ ■ 



+ Po 



where po, ■ ■ . , p r are polynomials and D represents the stan- 
dard derivation -4-. Such operators act in a natural way 
on elements of a differential ring containing the polynomi- 
als. An object y is called a solution of the operator if P 
applied to y yields zero. We are interested in finding the hy- 
perexponential solutions of a given operator. An object y is 
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called hyperexponential if the quotient D(y)/y can be identi- 
fied with a rational function. Typical examples are rational 
functions (e.g. (5x + 3) /(3a; + 5)), radicals (e.g. \/x + 1), ex- 
ponentials (e.g. exp(3x 2 —4) or exp(l/a;)), or combinations 
of these (e.g. \Jx + 1 exp(a; 9 / (x — 1))). Equivalently, y is 
called hyperexponential if there is some first order operator 
qiD + qo with qo,qi polynomials which maps y to zero, ff 
we regard differential operators as elements of an operator 
algebra C(x)[D\, then there is a one-to-one correspondence 
between the hyperexponential solutions y of an operator P 
and its first order right hand factors. In other words, if y is 
a hyperexponential term with (qiD + qo) ■ y — 0, then y is 
a solution of P if and only if there exist rational functions 
wo, . . . , u r _i such that 



P = 



.iZT- 1 + U T - 2 D r 



+ 



uo)(q\D + q ) 



Algorithms for finding the hyperexponential solutions of a 
linear differential equation (or equivalently, the first order 
right hand factors of the corresponding operators) are known 
since long. They are needed as subroutine in algorithms for 
factoring operators or for finding Liouvillean solutions. See 
Chapter 4 of [11] for details and references. 

Classical algorithms first compute "local solutions" at sin- 
gular points (cf. Section 2.3 below) and then test for each 
combination of local solutions whether it gives rise to a hy- 
perexponential solution. This leads to a combinatorial ex- 
plosion with exponential runtime. The situation is similar to 
classical algorithms for factoring polynomials over (Q, which 
first compute the irreducible factors modulo a prime and 
then test for each combination whether it gives rise to a 
factor in <Q[x]. 

The algorithm of van Hoeij [12] avoids the combinatorial 
explosion as follows. It picks one local solution and consid- 
ers the operator Q = qiD + qo with qi,qo £ C((x)) which 
annihilates it. This operator is a right factor of P, though 
not with rational coefficients. The algorithm then constructs 
(if possible) a left multiple B of Q with rational coefficients 
of order at most r — 1. This leads to a nontrivial factor- 
ization P = AB in C(x)[D]. The procedure is then ap- 
plied recursively to A and B until a complete factorization 
is found. The first order factors in this factorization give 
rise to at most r hyperexponential candidate solutions (pos- 
sibly up to multiplication by a rational function). These are 
then checked in a second step. Van Hoeij 's algorithm re- 
minds of the polynomial factorization algorithm of Lenstra, 
Lenstra, Lovasz [6, 15], which picks one modular factor and 
constructs (if possible) a multiple of this factor with integer 



coefficients but smaller degree than the original polynomial. 
This multiple is then a proper divisor in <Q[x]. 

The algorithm we propose below avoids the combinato- 
rial explosion in a different way. We start from the local 
solutions and regard them as asymptotic expansions of com- 
plex functions. By means of effective analytic continuation 
and arbitrary-precision numerical evaluation, we compute 
the values of these functions at some common ordinary ref- 
erence point. Then a linear algebra algorithm is used to 
determine a small list of possible combinations of local solu- 
tions that may give rise to hyperexponential ones, possibly 
up to multiplication by a rational function. These are then 
checked in a second step. Our approach was motivated by 
van Hoeij's polynomial factorization algorithm [14], which 
associates to every modular factor a certain vector and then 
uses lattice reduction to determine a small list of combina- 
tions that may give rise to proper factors. 

Although our algorithm avoids the combinatorial explo- 
sion problem, we do not claim that it runs in polynomial 
time. Indeed, no polynomial time algorithm can be expected 
because there are operators P which have hyperexponential 
solutions y that are exponentially larger than P. Also van 
Hoeij [12] makes no formal statement about the complexity 
of his algorithm. It is clear though that his algorithm is su- 
perior to the naive algorithm. Similarly, we believe that our 
algorithm has chances to outperform van Hoeij's algorithm, 
at least in examples that are not deliberately designed to ex- 
hibit worst case performance. The reason is partly that dur- 
ing the critical combination phase we only work with floating 
point numbers of moderate precision while van Hoeij's algo- 
rithm in general needs to do arithmetic in algebraic number 
fields whose degrees may grow during the computation. An- 
other advantage of our algorithm is that it is conceptually 
simpler than van Hoeij's, at least if we take for granted that 
we can compute high-precision evaluations of D-finite func- 
tions. 

2. PRELIMINARIES 

In this section, we recall some results from the literature 
and introduce notation that will be used in subsequent sec- 
tions. 

2.1 Differential Fields and Operator Algebras 

A differential ring/field is a pair (K,D) where if is a 
ring/field and D: K — > K is a derivation on K, i.e., a map 
satisfying D{a+b) = D(a)+D(b) and D(ab) = D(a)b+aD(b) 
for all a,b G K. Throughout this paper, we consider the dif- 
ferential field K = C(x), where C is some (computable) 
subfield of C, together with the derivation D : K — > K de- 
fined by D(c) = for all c G C and D(x) = 1. For simplicity, 
we assume throughout that C is algebraically closed. 

A differential ring/field E is called an extension of K if 
K C E, and the derivation of E restricted to K agrees with 
the derivation of K. 

By K\D] we denote the set of all polynomials in the in- 
determinate D with coefficients in K. Addition in K\D] is 
defined in the usual way, and multiplication is defined sub- 
ject to the commutation rule Da = aD + D(a) for a G K. 
The elements of K[D] are called operators, and they act 
on the elements of some extension E of K in the obvious 
way: If P — po + piD + ■ ■ ■ + p r D r is an operator of or- 
der r and y G E, then P ■ y := YH=iPi D ' l (v) G E - The 
noncommutative multiplication is compatible with operator 



application in the sense that we have (PQ) ■ y = P ■ (Q ■ y) 
for all P, Q G K[D] and all y G E. 

The elements y G E such that P ■ y = form a C-vector 
space V with dim V < r. By making E sufficiently large it 
can always be assumed that dim V = r. 

2.2 Hyperexponential Terms 

Let E be an extension of K. An element h G E \ {0} is 
called hyperexponential over K if D(h)/h G K. Equivalently, 
h is hyperexponential if Q- h = for some nonzero first order 
operator Q G JC[-D]- 

Two hyperexponential terms hi , /i2 are called equivalent 
if hxjhi G K. For example, the terms exp(3x 2 — a;) and 
(1 — 2a;) 2 exp(3a; 2 — x) are equivalent, but exp(3a; 2 — x) and 
(1 — 2x)^ exp(3x 2 — x) are not. (Here and below, we use 
standard calculus notation to refer to elements of some ex- 
tension E on which the derivation acts as the notation sug- 
gests, e.g. D(exp(3a; 2 — x)) = (6a; — 1) exp(3a; 2 — a;).) 

Every hyperexponential term can be written in the form 
h = exp(J v), where v is a rational function. The additive 
constant of the integral amounts to a multiplicative constant 
for h, which is irrelevant in our context, because P ■ h — 
if and only if P ■ (ch) = for every c G C \ {0}. If we con- 
sider the partial fraction decomposition of v and integrate 
it termwise, we obtain something of the form 

n 

9 + ^2li log (Pi) 
i=l 

with g G K, 71, . . . , 7„ G C and monic square free pairwise 
coprime polynomials pi G C[x]. In terms of this representa- 
tion, two hyperexponential terms are equivalent if the differ- 
ence of the corresponding rational functions g is a constant 
and any two corresponding coefficients 7; differ by an inte- 
ger. 

The equivalence class of a hyperexponential term h is 
called the exponential part of h. The motivation for this ter- 
minology is that when we are searching for some hyperexpo- 
nential solution h of P and we already know its equivalence 
class, then we can take an arbitrary element ho from this 
class and make an ansatz h — uho for some rational func- 
tion u G K. The operator P := P <g> (D - D [) / h h ° ) ) G K[D] 

then has the property that u is a solution of P if and only if 
uho is a solution of P. This reduces the problem to finding 
rational solutions, which is well understood and will not be 
discussed here [1, 11]. 

2.3 Local Solutions 

Consider an operator P G C(x)[D] of order r. By clear- 
ing denominators, if necessary, we may assume that P G 
C[x][D], say P = p r D r + ■■■ +po with p r 0. A point 
z G C U {00} is called singular if z is a root of p r , or z — 00. 
A point which is not singular is called ordinary. Note that 
there are only finitely many singular points, and that we 
include the "point at infinity" always among the singular 
points. 

If z = is an ordinary point then P admits r linearly inde- 
pendent power series solutions. If z = is a singular point, 
it is still possible to find r linearly independent generalized 
series solutions of the form 

m 

x a exp(u(x~ 1/s )) bk{x 1/s ) log(a;) fe (1) 

fe=0 



where a G C, u G C[x] with w(0) = 0, s G N, m G N 
and bo,..., b m G C[[x]]. We call these solutions the local 
solutions at 0. The computation of such solutions is well- 
known and will not be discussed here [13, 11]. 

Two series as in (1) are called equivalent if they have the 
same u and s and the difference of the respective values of 
a is in The equivalence classes of generalized series un- 
der this equivalence relation are called the exponential parts 
of the series. Adopting van Hoeij's notation and defining 
Exp(e) := exp(J ^) for e G C[s _1,/s ], we have that Exp(ei) 
and Exp(e2) are equivalent iff ei — e2 G \7L. Note that if 
m — and s = 1, two series are equivalent iff their quotient 
can be identified with a formal Laurent series. We will from 
now on make no notational distinction between Exp(e) and 
its equivalence class. 

A point 2^0 can be moved to the origin by the change 
of variables x = x — z (if z G C) or x = 1/x (if z = oo). If 
P is the operator obtained from P by replacing x by x + z 
or 1/x, then a local solution of P G C[a:][D] at z is defined 
as the local solution of P G C[x][ZJ] at 0. 

Throughout the rest of this paper, we will use the follow- 
ing notation. P is some operator in C[x][_D] of order r, by 
zi, . . . , z n -i G C we denote its finite singular points, z n = oo. 
We write Xi = x — Zi (i = 1, . . . , n — 1) and x n = 1/x for 
the variables with respect to which the singularities at Zi ap- 
pear at the origin. For i = 1, . . . , n, we consider the vector 
space Vi generated by all local solutions at Zi. There may 
be solutions with different exponential parts, say ii different 

— 1/S' ' 

parts Exp(ej,i), . . . ,Exp(ei,< i ) for e»,j G C\x i ,,3 \. By 

Vi,j = Vi nExp(e lJ )C((^ 1/Sl ' 3 ))[loga; ! ] 

we denote the vector space of all local solutions of P at Zi 
with exponential part (equivalent to) Exp(e ii £ i ). Our Vi,j 
are written V ei } (P) in van Hoeij's papers [13, 12]. 

The condition in the definition of equivalence that the 
difference of corresponding values of a be an integer (rather 
than, say, requiring exactly the same value of a) ensures 
that the Vij are indeed vector spaces, because if some Vij 
contains, for example, the two series 

x a (l + x + x 2 + ■ ■ • ) and x a (l + x + 3x 2 + ■ ■ ■ ) 

then it must also contain their difference x a (2x 2 + ■ ■ ■ ) = 
x <*+ 2 (2 + ■■■). 

2.4 Analytic Solutions 

It is classical that the formal power series solutions y of P 
at an ordinary point z G C actually converge in a neighbour- 
hood of z and thus give rise to analytic function solutions y 
of P. The correspondence is one-to-one. For any other or- 
dinary point z' G C and a path z ~» z' avoiding singular 
points there exists a matrix M z ^, z * G <D rxr such that 

(& y (z')y-j o = M^ z ,(D> V (z)y 3 :i 

for every solution y analytic near z. There are algorithms [4, 
8] for efficiently computing the entries of M z ^, z i for any 
given polygon path z "~> z' with vertices in <Q to any de- 
sired precision. In other words, we can compute arbitrary 
precision approximations of y and its derivatives at every 
ordinary point ("effective analytic continuation"). 

Assume now that is a singular point, and consider the 
case s = 1 and m = 0, i.e., let y — Exp(e)6 for some e G 
C[a; _1 ] and 6 G C[[x]] be a formal solution of P. To give an 



analytic meaning to Exp(e) = exp(J f ) = exp(«+a log x) — 
x a exp(it) (for suitable a G C and u G C[a; _1 ]) amounts to 
making a choice for a branch of the logarithm. Every choice 
gives rise to the same function up to some multiplicative 
constant. 

Since Exp(e)6 is a solution of P iff b is a solution of the 
operator P <g> (D + -), we may assume that e = 0. Then 
the problem remains that the formal power series y — b 
may not be convergent if is a singular point. However, by 
resummation theory [2, 3] it is still possible to associate to 
y an analytic function y defined on some sector 

A = A(d, ip, p) := {z e € : < \z\ < p A \d - arg«| < ip/2} 

(with d G [0, 2n], p,cp > 0) such that y is the asymptotic 
expansion of y for z — > in A. 

The precise formulation of this result is technical and not 
really needed for our purpose (see [3, Chap. 6, 10, and 11] or 
[2, Chap. 5-7] for full details). It will be more than sufficient 
to know the following facts: 

• For every k — (ki, . . . , k q ) G with ki > •• • > k q 
and every d = (di, . . . , d q ) G [0, 27r] 9 such that 

\dj+i -dj\< (fc~ + \ - fc" 1 )! , j = 1, • • • ,q - 1, 

one constructs [3, §10.2] a differential subring <C{x}k,d 
of C[[x]] [3, Theorems 51 and 53] which contains the 
ring C{a;} of all convergent power series. 

• There is a differential ring homomorphism [3, The- 
orems 51 and 53] Sk,d from <C{x}k,d to the germs 
of analytic functions defined on sectors of the form 
A(d\, <fi, p) for suitable tp, p > 0, with the property that 
for every y G < C{x}k,d the function Sk,d{y) has y as its 
asymptotic expansion for z — > [3, §10.2, Exercice 2]. 
The Sk,d map convergent formal power series to their 
sum in the usual sense [3, Lemmas 8 and 20]. 

• For a given operator P G C[x][D] of order r, one 
can compute a tuple fc and finite subsets T>i , . . . , T> q 
of [0, 27r] such that any y G C[[x]] with P-y = belongs 
to C{x}k,d for all d as above with d\ (/ T>\, . . . , d q £ 
T> q . Additionally, given such a d, one can compute 
if, p > such that each Sk,d(y) is defined on A(di, tp, p). 

• Furthermore, given a point z G A(di,y>, p), a precision 
e > 0, and y G C[[x]] with P-y = 0, one can efficiently 
compute an approximation Y e of the vector Y(z) — 
{D j S k ,d{y)) r 3 Zl such that \\Y(z) -Y e \\< e. 

The computational part of the last two items is a special 
case of Theorem 7 of van der Hoeven [10]. As an application, 
van der Hoeven [9] shows how to factor differential operators 
using numerical evaluation. Note that our kj correspond to 
1/kj in van der Hoeven's articles, and the components of the 
tuples k and d appear in reverse order. 

Also observe that in the last item, z is an ordinary point, 
so that from there we can use effective analytic continuation 
to compute values of Sk,d(y) and its derivatives at any other 
ordinary point. 

3. OUTLINE OF THE ALGORITHM 

A hyperexponential term h can be expanded as a gener- 
alized series at every point z G C U {oo}, in particular at 
its singularities. The resulting generalized series are local 



solutions of P if h is a solution of P. If h — exp(J v) is a 
hyperexponential solution where v £ C(x), and if we write 
the partial fraction decomposition of v in the form 

ei e 2 e n 
u = 1 h h 77-, 

a; — zi a; — 2:2 1/a; 

where the e< are polynomials in x" 1 , then expanding this h 
at «j yields a generalized series in Xi whose exponential part 
matches Exp(ej). The components e» in the decomposition 
of v must hence show up among the exponential parts of the 
local solutions of P. 

If Exp(ej,i), . . . , Exp(ej i £ i ) are (representatives of) the dif- 
ferent exponential parts that appear among the local so- 
lutions at Zi, then any hyperexponential solution must be 
equivalent to the term exp(J( + ■ ■ ■ + L '~' 3 " )) for some 
tuple (ji, . . . ,j„). It then remains to check for each of these 
candidates whether some element of its equivalence class 
solves the given equation. The basic structure of the algo- 
rithm for finding hyperexponential solutions is thus as fol- 
lows. 

Algorithm 1. Input: a linear differential operator P = 
po + piD + ■ ■ ■ + p,D r , p r / 0, with coefficients in C[x]. 
Output: all the hyperexponential terms h with P ■ h = 0. 

1. Let zi,. ■ ■ , z n _i £ C be the roots of p r in <D, and let 
z n =co. 

2. For i = 1, . . . , n do 

3. Find the exponential parts Exp(ei,i), . . . , Exp(ei/ i ) 
of the local solutions of P at z%. 

4- Determine a set U C {1,...,£\} x ■■■ x {1,...,£„} 
s.t. for every hyperexponential solution h equivalent to 

ex P(J Z)t=i %f") we have 0'i> ■ • • >■?») e u - 

5. For each (ji, . . . , j n ) £ U do 

6. Let ho := ex p(/X/"-i ant ^ com P u ^ e ^ e QJ?er- 
ator P:=P®(D- D( ^ ( f ) ■ 

7. Compute a basis {ui, . . . , u m } C C(x) o/ i/ie vector 
space of all rational solutions of P, and output uiHq, 
. . . , u m ho. 

There is some freedom in step 4 of this algorithm. A naive 
approach would simply be to take all possible combinations, 
i.e., U = {l,...,£i} x ■•■ x {l,...,£ n }- This is a finite 
set, but its size is in general exponential in the number of 
singular points. For finding a smaller set U, Cluzeau and 
van Hoeij [5] use modular techniques to quickly discard un- 
necessary tuples. Our algorithm, explained in the following 
section, addresses the same issue. It computes a set U of at 
most r tuples. 

4. THE COMBINATION PHASE 

In general, the differential operator P may have several 
different solutions with the same exponential part, i.e., the 
dimension of the vector spaces Vij might be greater than 
one. In this case, it might be that Vij contains some series 
which is the expansion of a hyperexponential solution h at Zi 
as well as some other series which are not. If we compute 
some basis of Vij, we cannot expect it to contain the expan- 
sion of h. Instead, each basis element will in general be the 
linear combination of this series and some other one. Now, if 
the expansion of h at some other singular point z^ belongs 



to the space Vi\j' (which possibly also has higher dimen- 
sion), then, in some sense, h must belong to the intersection 
of the vector spaces V^j and Vi/ji. 

Our algorithm is based on testing which intersections are 
nontrivial. To make these intersections meaningful, we must 
first map the vector spaces we want to intersect into a com- 
mon ambient space W . Let E be some differential ring con- 
taining C(x) as well as all the hyperexponential solutions 
of P, and let W C E be the C-vector space generated by 
solutions of P in E. For each i, let Hi be some vector space 
homomorphism 

It 

0Exp( eiJ )C((^ 1/Sl ' J ))[log^] DVi^UW 

3 = 1 

with the following properties: 

1. The sum 7Ti(Vj,i) + ■ ■ ■ + ^(V^) is direct. 

2. If h £ W is hyperexponential, then n~ (h) contains 
the formal series expansion h of h at Zi, possibly up to 
a multiplicative constant. 

Define Wij := ffi(Vi : j). If h is some hyperexponential solu- 
tion of P, say with exponential part 

\J \ Xl x 2 x n J J 

then h £ Wi,j t for all i, and hence the vector space Wi t j 1 n 
■ ■ ■ H Wn,j n is not the zero subspace (because it contains at 
least h). Our main observation is that there can be at most 
r tuples j = (ji, ... ,j„) for which Wj ^ {0}, and that they 
can be computed efficiently once we have bases of the Wij. 

Postponing the discussion of making the Hi constructive 
to the next section, assume for the moment that W is some 
vector space over C, let r = dim W < oo be its dimension, 
and suppose we are given n different decompositions of sub- 
spaces of W into direct sums: 

e Wi# e • • • e w Ml c w, 
w-2,1 e w 2 ,2 e ■ • ■ e w 2 ,t 2 c w, 

w n ,i e w„, 2 e • • • e w n ,i n c w. 

Without loss of generality, we may make the following as- 
sumptions: 

• Each direct sum ^ Wij is in fact equal to W. If 
not, add one more vector space to the sum. 

• £\ — £2 — ■ ■ ■ = In =: £. If not, pad the sum with 
several copies of {0}. 

• £ < r. If not, then because the sums are supposed 
to be direct, each decomposition must contain at least 
£ — r copies of {0}, which can be dropped. 

Lemma 2. There are at most dim W = r different tuples 

i = (ji,...,in)£{i,...,n n 

such that Wj := Wi, jl n W 2 j 2 H •• • H W„,j„ {0}. 

Proof. Induction on n. For n = 1, there are only £ < r 
different tuples altogether: (1), (2), ... , (£), so the claim is 



obviously true. Suppose now that the claim is shown for 
the case when n — 1 decompositions of some vector space 
are given. Let U C {1, ... ,£} n be a set of tuples j with 
Wj 7^ {0}. Partition the elements of U according to their 
first components, 

U = Ui U U 2 U ■■• U U e , 

i.e., Uk is the set of all tuples j whose first component is k, 
for k = 1,...,£. 

For all j = (k,j 2 , ■ ■ ■ , j n ) G U k we have {0} Wj C Wi,*,. 
Therefore, (j2, ■ ■ ■ , jn) G {1, . . . ,£} n_1 is a valid solution tu- 
ple for the modified problem with W(j := Wi+ij PI 
(i = 1, . . . , n — 1, j = !,...,£) in place of (i = 1, . . . , n, 
j — 1, . . . ,£). By induction hypothesis, since the W(j form 
n — 1 decompositions of the space Wi,*,, there are at most 
dimWi ifc tuples {j 2 ,...,j n ) with Wfa,...,j m ) {0}- Con- 
sequently, there are altogether at most dimWi.fc = 
dim W — r different tuples for the original space W . m 

The desired index tuples can be computed efficiently using 
dynamic programming, as shown in the following algorithm. 

Algorithm 3. Input: a vector space W of dimension r, and 
a collection of subspaces Wi,j (i = l,...,n; j = 1, ...,£) 
such that W = © J= i Wijj f or ' = 1, • - • , n an ^ £ <r. 
Output: the set U of all tuples j = (ji, ■ ■ ■ ,jn) with the 
property Wj = f|™ =1 Wi, it + {0}. 

1. U:={ (j) : Wi,j j4 {0} } 

2. For z = 2, . . . , n rfo 

5- Ufi ew • — 

4. For j = !,...,£ do 

5. For k e U do 

6. If W k n Wi,j ± {0} then 

7. U new := Unao U {append(fe, j)} 

8. U . — Unew 

9. Return U 

Theorem 4. Algorithm 3 is correct and needs no more than 
8nr 4 operations in C , if the bases of the Wk are cached. 

Proof. Correctness is obvious by line 6 and the fact that 
whenever k — (fei, . . . , k n ) is such that Wk 7^ {0} then we 
necessarily also have W(k!,...,k _-,) 7^ {0}. 

For the complexity, we first show that it is a loop invari- 
ant that Wk t l~l Wk 2 = {0} for any two distinct k\,k 2 G U. 
This is clear for i = 1 by line 1 and the assumption in 
the algorithm specification that W = Wi,j is a direct 

sum. Assume it is true for some i and consider the situa- 
tion right before line 8. At this point, for any two distinct 
tuples k\,k 2 G U we have Wk x n Wk 2 = {0} by induction 
hypothesis. We have to show that the same is true for any 
two distinct tuples ki,k 2 G U n ew By line 7, any such tu- 
ples have the form fei = (u\,ji), k 2 = (u 2 ,j 2 ) for some 
ui,u 2 G U and ji,j 2 G {1, ...,£}. The tuples ki,k 2 are 
distinct if ui 7^ u 2 or ji 7^ j 2 . If ui 7^ u 2 , then by induction 
hypothesis W U1 n W 7 ^ = {0}, and therefore also 

w kl n w k2 = (w*! n w i:jl ) n (w„ 2 n w t , j2 ) 
= {0} n n w id2 = {0}. 



Similarly, if ji / j'2, then Wij ± CiWij 2 = {0} by the assump- 
tion that W — Wij is a direct sum. Therefore 

w kl n w k2 = (w ai n n (W U2 n w 4 , ia ) 
= W Ul n w„ 2 n {0} = {0}. 

This completes the proof of the loop invariant fei 7^ k 2 =>■ 

VF fel n w k2 = {o}. 

A consequence of this invariant is that ^2 keu dim Wk < r 
in every iteration. Since the sum W = W»j is direct, 

we also have •_ 1 dim VFi, j < r in every iteration. The 
intersection of two subspaces of W of dimensions di , d 2 can 
be computed using no more than 

min(r, di + d 2 ) 2 max(r, di + d 2 ) 

operations in C. For the total cost of the algorithm we 
therefore obtain, writing Ui for the set U in the ith iteration 

dim Wi, 

^ ^ min(r, d k + dij) 2 

i=2 3=1 fee(7; 



i=2 j=l 



< 2r^2(£r 2 + 2r 3 + r 3 ) 

< 8nr 4 



i=2 
4 



In the second step, we have used the bounds ^2 keu . d\ < 
r 2 and \Ui\ < r, which follow from ^2 keu .dk < r and 
Lemma 2, respectively. In the third step, we used the bound 
X^-i d%,j < t ' 2 j which follows from X^=i — r - ■ 

If the objective is just to show that the algorithm runs in 
polynomial time, a simpler argument applies. It suffices to 
observe that all the intersections can be done with a number 
of operations which is at most cubic in r, then taking also 
into account that we always have \ U\ < r by Lemma 2, the 
bound 0(n£r 4 ) — 0(nr 5 ) follows immediately. 

5. NUMERICAL EVALUATION AT A 
REFERENCE POINT 

We now turn to the question of how to construct the mor- 
phisms -Ki. The basic idea is to choose a reference point zo 
that is an ordinary point of P, and let W be the space of 
analytic solutions of the equation in a neighborhood of zo. 

For each singular point z%, let Aj be a 
sector rooted at Zj for which all formal 
power series appearing in the generalized 
series solutions of P at Zi admit an inter- 
pretation as analytic functions via some 
operator Sk,d (depending on i, but not 
on the series), as described in Section 2.4. 
Such sectors exist and can be computed 
explicitly. Next, let ji (i — l,...,n) be polygonal paths 




from Zi to zq avoiding singular points and leaving the start- 
point through Aj (meaning that for some e > all the points 
on ji with a distance to Zi less than e should belong to Aj). 
Such paths exist. The analytic interpretations of the gener- 
alized series solutions at the singular points z; defined in A; 
admit a unique analytic continuation along the paths ji to 
the neighborhood of zq. 

We define Tv t : Vj — > W as follows. Let Vf j be the subspace 
of Vij consisting of generalized series (1) with s = 1 and m = 
0, and let V{ a be a linear complement of Va°a in Vij. If y G 
V°a i.e., if y = Exp(e ij )fe with e itj G C^ 1 ] and b G C[[xi]], 
define TTi(y) to be the unique analytic continuation of the 
function E(ej,j)>Sfe,d(6) along ji to zo, where E(eij) refers to 
the function z h-> exp(J z eij/xi) with some arbitrary but 
fixed choice of the branch of the logarithm, and Sk,d is as 
described in Section 2.4. Set iVi(y) = for y G V^'j, and then 
extend 7Ti to K by linearity. The precise values of Hi(Vi : j) 
depend on the choice of A, and d (which is arbitrary, within 
the limits indicated in Section 2.4), but, as shown below, the 
properties of these spaces used in the algorithm do not. 

Proposition 5. The functions Hi defined above satisfy the 
two requirements imposed in Section 4-' (1) i"<(Vi,i) + • • • + 
Kiiyi,^) is a direct sum; (2) if ft is a hyperexponential term, 
then 7r~ 1 (/i) contains the formal series expansion of h at zi, 
possibly up to a multiplicative constant. 

Proof. 1. Without loss of generality, we assume z; = 0. Let 
Vj G Vi,j (j = l,...,£i) and consider y = J2/=i Hi- Write 
ijj = x aj exp(uj)bj + y'j where y'j G Vl a, the (ctj,Uj) are 
pairwise distinct, Uj(0) = 0, and bj(0) 7^ unless the series 
bj is zero. Writing u~ = ^2 k u j,kX~ k , choose a direction 6 

such that pe l8 G Aj for small p and any two unequal Uj£e l9 
have different real parts. 

By changing x to e~ lB x, we can assume that d = 0. This 
tranforms Uj into ^ fc (u :)i fce lfce )x~ fc , so that the real parts of 
two polynomials Uj can be the same only if the Uj themselves 
are equal. Hence, we can reorder the nonzero terms in the 
expression of y by asymptotic growth rate, in such a way that 
the nonzero terms come first, ui = • • • = Ut and Recii = 
• ■ ■ = Re at , while 

z i e lw 2>z p e pw , z — > 0, z > 

for all p > t + 1 such that y p 7^ 0. Using the definition of -Ki 
and the fact that Sk,d{bj)(z) tends to bj(0) as z -> in the 
positive reals, it follows that 

t 

z- Rcai cxp(- Ml (z))^(y)(z) =^cA'(0K ImQy +o(l) 

i=l 

(as z — > 0, z > 0) for some nonzero constants Cj. Since 
the (aj , Uj ) are pairwise distinct by assumption and the 
(Keaj,Uj) are equal for j — 1, . . . ,t, the Im«j are pairwise 
distinct for j = 1, . . . , t. 

Now assume that 7Ti(y) = 0. Then, for all A > 0, the 
expression X^=i c 3^3'(0)(^ 2 ) mQ! 3 tends to as z — > 0, z > 
0. Choosing A = e p for p = 1, . . . , t, it follows that if not all 
the bj(0) were zero, the t x t determinant 

det((e p z) lImQ ') p q = z iIra( "' "' ; det((e lIma ') p ) j3 

would tend to zero as well, which however is not the case. 
Therefore bj (0) = for j = 1, . . . , t, and therefore ijj = for 
j — 1, . . . ,t, and therefore j/j = for j = 1, . . . , ii. 



2. Let h G W be hyperexponential. Then the expansion 
/i of ft at Zi is clearly a local solution, so ft 6 Vi,j for some j. 
We show that 7Ti(ft) = eft for some c G C The map 7Tj is 
a differential homomorphism because Sk,d is (as remarked 
in Section 2.4) and the (formal) exponential parts Exp(ej,j) 
are mapped to analytic functions satisfying the same differ- 
ential equations. Since ft is hyperexponential, it satisfies a 
first order linear differential equation. Since ft is the expan- 
sion of ft, it satisfies the same equations as ft. Since ni is 
a differential homomorphism, in(h) satisfies the same equa- 
tions as ft. Hence Hi(h) and ft satisfy the same first-order 
differential equation. The claim follows. ■ 

The definition of the maps 7Ti , . . . , 7r m as outlined above 
relies on analytic continuation, a concept which is only avail- 
able if C = (D. For actual computations, we must work in a 
computable coefficient domain. At this point, we use numer- 
ical approximations. By van der Hoeven's result quoted in 
Section 2.4, we are able to compute for every given y G Vi.j 
and every given e > a vector Y E G <Q(i) r with 

||K^)(-o))^-nL< £ . 

Using these approximations, the linear algebra parts of Al- 
gorithm 3 are then performed with ball arithmetic to keep 
track of accumulating errors during the calculations. The 
test in line 6 of this algorithm requires to check whether a 
certain matrix has full rank. There are two possible out- 
comes: If during the Gaussian elimination we can find in 
every iteration an entry which is definitely different from 
zero, then the rank of the matrix is definitely maximal and 
the intersection of the vector spaces is definitely empty. We 
are then entitled to discard the possible extension of the par- 
tial tuple under consideration. On the other hand, if during 
the Gaussian elimination we encounter a column in which all 
the entries are balls that contain zero, this can either mean 
that the intersection is really nonempty, or that the accuracy 
of the approximation was insufficient. In this case, in order 
to be on the safe side, we must consider the intersection as 
nonempty and include the corresponding tuple. 

Regardless of which initial accuracy e is used, this variant 
of Algorithm 3 produces a set of tuples that is guaranteed to 
contain all correct ones, but may possibly contain additional 
ones. With sufficiently high precision, the number of tuples 
in the output that actually have an empty intersection will 
drop to zero. We don't need to know in advance which 
precision is sufficient in this sense, because it is not dramatic 
to have some extra tuples in the output as long as they are 
not too many. As a pragmatic strategy balancing precision 
and output size, one might start the algorithm with some 
fixed precision e and let it abort and restart with doubled 
precision whenever \U\ exceeds 2r, say. 

Observe that the numerical approximation is only used to 
determine the tuple set U, and we do not use it to somehow 
reconstruct the exact symbolic hyperexponential solutions 
from it. We therefore don't expect to need very high preci- 
sion in typical situations. 

6. A DETAILED EXAMPLE 

Consider the operator 

P = Po+PiD + p 2 D 2 +p 3 D 3 G Q[x][D] 

where 



-is 



p„ = -105a; 20 + 3570a: 19 - 58026a; 18 + 556216a; 17 - 3456830a; 16 + 
14810744a; 15 - 45667732a; 14 + 104614932a; 13 - 182764261a; 12 + 
249940430a 11 - 276371642a; 10 + 257839924a- 9 - 211785148a; 8 + 
154714472a; 7 - 95675216a 6 + 45214304a 5 - 13863936a; 4 + 
1685888a 3 + 424960a 2 - 182784a; + 20480, 

16 



' + 



P1 = (x - l)a(105a; 19 - 3150a; 18 + 51456a 17 - 489796a; 16 
2938210a; 15 - 11903624a 14 + 34247824a 13 - 72603516a; 12 
116974957a 11 - 148046826a 10 + 153582952a; 9 - 137261696a; + 
109046080a; 7 - 75250624a 6 + 41559168a 5 - 16084864a; 4 + 
3278080a 3 + 163840a 2 - 231424a + 32768), 

p 2 = -i(x - 2) 2 (a - l) 3 a; 2 (30a; 15 - 693a; 14 + 7314a; 13 - 42905a 12 + 

155930X 11 - 378483a 10 + 649718a 9 - 828795a; 8 + 820160a; 7 - 
645092a; 6 + 398200a 5 - 182384a 4 + 54656a 3 - 5696a; 2 - 2944a; + 1024), 

p 3 = 4(a - 2) 4 (a; - l) 5 a 4 (15a 10 - 258a 9 + 1492a; 8 - 4446a 7 + 
8309a 6 - 10972a; 5 + 10520a 4 - 6456a 3 + 1552a 2 + 480a; - 256). 



The leading coefficient P3 has 13 distinct roots in (D, but 
those coming from the degree- 10-factor turn out to be ap- 
parent, so we can ignore them, ft thus remains to study the 
singular points z\ := 0, Z2 := 1, zz ■= 2, and 2:4 := 00. 

For each singular point, we find three linearly indepen- 
dent generalized series solutions with two distinct exponen- 
tial parts: 



Vi.i 


= Cfc,i 


Vl,2 


= Cyi, 2 


+ Cyi.3 


1/2,1 


= C02.1 


V 2 ,2 


= Cj/2,2 


+ Cj/2,3 


V 3 ,i 


= Cy 3 ,i 


V 3 ,2 


= Cj/3,2 


+ Cj/3,3 


Va,i 


= Cj/4,1 


1/4,2 


= Cl/4,2 


+ Cj/4,3 



where 
2/i,i 
2/1,2 

#1,3 



exp 



\fx(l — 5 



± x + Sl x 2 . 

g X T 32 X 



24 ~ 



384 ^ ~ 



f- ( 2 _ 7 3 , 9 4 
"\J X I X ^ I 22 



#2,1 = (a; - l) 3 + (x - l) 5 - |(x - l) 6 + ■ 



3/2,2 = exp 
j/2,3 = exp 

273,1 

2/3,2 
2/3,3 

2/4,1 
2/4,2 
2/4,2 



+ ■ 



(^)(l +§(*-!) + 

(_i r) ( (x _l)=» + i(a! _l ) 3 + ...) ) 

.l-f(s-2) + §(z-2) 2 -|i(*-2) 3 + . 

F ^ F exp(^)(l + ii(x-2)+. 
: F ^ F exp(^)((x-2) 3 + I( a; -2) 4 + . 

■■ x(l + 3a;" 1 + 9a;" 2 + ^a;" 3 + 74a;" 4 + ■ • 
■■y/x(l + x- 1 + fa;" 2 + fx" 3 + ■ 
y/x(x 3 + x + fx' 1 + 4jf a;" 2 + ■ 



Let us choose zo = 3 as ordinary reference point and take 
the branch of the logarithm for which <Jx is positive and 
real on the positive real axis. The example was chosen in 
such a way that all the power series are convergent in some 
neighborhood of the expansion point, so that we do not need 
to worry about sectors and resummation theory but can 
use the somewhat simpler algorithm for effective analytic 
continuation in the ordinary case to compute the values of 
the analytic functions y%j := TTi(i)i,j) {i = 1,...,4; j — 



1,2,3). The vectors (yi,j(zo), Dyij(zo), D 2 yi,j(z )) to five 
decimal digits of accuracy are as follows. 



Wi, 



W 3 , 



Wa 



'-200.15N 

322.46 
,-1184.8, 

' 30.349 N 
-48.896 
v 179.66 j 

' .74285 \ 
-.061904 
v .14960 / 

' 30.349 N 
-48.896 
> 179.66 , 



W 2 ,2 



W 4 ,: 



/ -70.513\ 
-46.308 . 
^-101.17/ 

'12.494N , 
5.2891 , 
v 13.066/ ' 

' 15.580 \ 
-31.307 . 
v, 105.26 / 

' 2.8557 \ 
-.23797 . 
v .57510 / 



-156.55^ 

-91.322 

-205.47/ 



'77.105^* 
44.216 
v 99.931 j 

/4.5433X 
, 2.6503 
\5.9631/ 

/63.199X 
, 41.308 
V90.353/ 



We now go through Algorithm 3. Start with the partial 
tuples (1) and (2) corresponding to the vector spaces W\ t \ 
and Wi t 2, respectively. To compute the intersection of Wi,\ 
and W-j,! we apply Gaussian elimination to the 3 x 2-matrix 
whose columns are the generators of Wi,i and Wa.i: 



-200.15 
322.46 
-1184.8 



30.349 N 
-48.896 
179.66 i 



-200.15 30.349^ 
0.00 
0.00 / 



The notation 0.00 refers to some complex number z with 
\z\ < 5-10" 3 , which may or may not be zero, while the blank 
entries in the left column signify exact zeros that have been 
produced by the elimination. As the remaining submatrix 
does not contain any entry which is certainly nonzero, we re- 
gard the intersection as nonempty, which in this case means 
Wi,i = W 2 ,i- The partial tuple (1) is extended to (1,1). 

The intersections Wi,i PI W 2 ,2 and Wi,2 PI W2,i turn out 
to be trivial, as they have to be if we really have Wi,i = 
Wi,2, because the sums Wi,i © Wi,2 and Wz,i © Wa,2 are 
direct. It thus remains to consider the intersection Wi,2 H 
W2,2. Applying Gaussian elimination to the 3 x 4-matrix 
whose columns are the generators of Wi,2 and W2,2, we find 

/-70.513 -25.596 12.494 77.105X 

-46.308 2.1330 5.2891 44.216 
V-101.17 -5.1548 13.066 99.931/ 

/-70.513 -25.596 12.494 77.105X 
— ► -17.50 4.440 9.777 , 

V 0.00 0.00 / 

which suggests that we have W\,2 = W2.2- We extend the 
partial tuple (2) to (2,2). At the end of the first iteration, 
we have U = {(1, 1), (2, 2)}. 

In the second iteration, we find W(i i) PI W3,i = {0} and 
W/(i,i) Q Wa,2, so we extend the partial tuple (1,1) to 
(1,1,2) and record H/(i,i,2) = W(i,i) = Wi,i. Furthermore 
we find Wz,\ C W(2,2), so we extend (2,2) to (2,2,1) and 
record W[2,2,\) = W^.i- Finally, there is a nontrivial inter- 
section between W( 2y 2) and W3,2- 

/12.494 77.105 15.580 4.5433X 

5.2891 44.216 -31.307 2.6503 
\13.066 99.931 105.26 5.9631/ 

/12.494 77.105 15.580 4.5433X 
— > -27.34 89.53 -1.72 

V 216. 0.00 / 



suggests a common subspace of dimension 1 generated by 
the second listed generator of Wz,2- We therefore extend 
the partial tuple (2, 2) to (2, 2, 2) and record Wt2fi,2) — 
[(4.5433, 2.6503, 5.9631)]. At the end of the second iteration, 
we have U = {(1,1, 2), (2, 2,1), (2, 2, 2)}. 

For the final iteration, we see by inspection that W4,i = 
Wi t \ — W^ni^), so we extend (1, 1, 2) to (1, 1, 2, 1). Because 
dimVK4,i = 1 and the sums of the vector spaces are direct, 
the other two partial tuples cannot also have a nontrivial 



intersection with W4,i, nor can W t 



ial. 



(1,1,2) 



n Wi 2 be nontriv- 



We do however have Wt 2,2,1) 



W4,2 and W( 2 ,2,2) ^ 
W4 $ 2, so the algorithm terminates with the output U = 
{(1,1, 2,1), (2, 2, 1,2), (2, 2, 2, 2)}. 

At this point we know that every hyperexponential solu- 
tion of the operator P must have one of the following three 
exponential parts: 



__!_ex ( l \ 1 

(a; — 2) 2 ^ Va; x — 



^ exp (^l) 



V^exp 



x - 1 x-2 



from (1,1,2,1) 
from (2,2,1,2) 
from (2,2,2,2). 



Following the steps of Algorithm 1, it remains to check 
whether some rational function multiples of these terms are 
solutions of P. The important point is that we have to do 
this only for three different candidates, while the naive al- 
gorithm would have to go through all 2 4 = 16 combinations. 
Indeed, it turns out that P has the following three hyperex- 
ponential solutions: 



{x-2) 2 

(x — 2)x 2 \fx expf - 

V x — 1 



(x-\f (\ 1 \ r- I 

■ exp I — I -I, y/xexp - 

Va; x — 2J \: 



+ 



x-2 



7. CONCLUDING REMARKS 

Our algorithm as described above takes advantage of the 
fact that series expansions of hyperexponential terms can- 
not involve exponential terms with ramification (s > 1) or 
logarithms (m > 0), by letting the morphisms m map all 
these irrelevant series solutions to zero. As a result, we get 
smaller vector spaces Wij, which not only reduces the ex- 
pected computation time per vector space intersection but 
also makes it somehow more likely for intersections to be 
empty, thus decreasing the chances of getting tuples that do 
not correspond to hyperexponential solutions. 

As a further refinement in this direction, it would be desir- 
able to exploit the fact that if h — Exp(e)& is the expansion 
of some hyperexponential term h, then the formal power 
series b must be convergent in some neighborhood of the ex- 
pansion point. Instead of the vector spaces Wij used above, 
it would be sufficient to consider the subspaces W(j C Wij 
corresponding to generalized series solutions involving only 
convergent power series. Besides the advantage of having 
to work with even smaller vector spaces, an additional ad- 
vantage would be that the numerical evaluation becomes 
simpler because algorithms for the regular case [4, 8] be- 
come applicable. Implementations of these algorithms are 
available [7], which to our knowledge is not yet the case for 
van der Hoeven's general algorithm for the divergent case 
[10]. Unfortunately however, it is not obvious how to com- 
pute from a given basis of Wi,j a basis of the subspace W(j. 



Miller's algorithm [16] numerically solves a similar problem, 
but so far we have not been able to turn the underlying con- 
vergence statements into explicit error bounds that would 
yield an algorithm producing output with certified precision. 

Finally, it would of course be also interesting to see an 
analog of our algorithm for finding hypergeometric solutions 
of linear recurrence equations with polynomial coefficients. 
A translation is not immediate because there is no notion of 
local solution around a finite singularity in this case. 



8. 
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